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Abstract 

Due to the sign problem, it is exponentially difficult to study QCD on the lattice at finite chemical potential. We 
propose a method -an overlap improving multi-parameter reweighting technique- to alleviate this problem. We apply 
this method and give the phase diagram of four-flavor QCD obtained on lattices 4 4 and 4 ■ 6 3 . Our results are based 
on O(10 3 - 10 4 ) configurations. 



1 Introduction 



Quantum Chromodynamics (QCD) at finite temperature (T) and/or chemical potential (fj,) is of fundamental importance, 
C^) . since it describes relevant features of particle physics in the early universe, in neutron stars and in heavy ion collisions. 
' According to the standard picture, as baryon density rises there is a change from a state dominated by hadrons -protons 
and neutrons- to a state dominated by partons -quarks and gluons. In addition, recently a particularly interesting, rich 
phase structure has been conjectured for QCD at finite T and \i [0, |^]. Of immediate interest is the existence and the 
location of the critical point in the T-fi plane in three-flavor QCD, since it can be explored by heavy-ion experiments. 
Clearly, the transition from hadronic to partonic (or to superconducting/superfluid) state is a fully non-perturbative 
phenomenon. 

Lattice gauge theory is a reliable systematic technique to study the non-perturbative features of QCD. QCD at finite 
fj, can be formulated on the lattice ||; however, standard Monte-Carlo techniques can not be used at (J, ^ 0. The reason is 
that for non- vanishing real fi the functional measure thus, the determinant of the Euclidean Dirac operator- is complex. 
This fact spoils any Monte-Carlo technique based on importance sampling. 
j^J ' Several suggestions were studied in detail to solve the problem. We list a few of them. 

In the large gauge coupling limit a monomer-dimer algorithm was used For small gauge coupling an attractive 
' approach is the "Glasgow method" ^ in which the partition function is expanded in powers of exp(/i/T) by using 
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an ensemble of configurations weighted by the fj,=0 action. After collecting more than 20 million configurations only 
unphysical results were obtained: a premature onset transition. The reason is that the /i=0 ensemble does not overlap 
sufficiently with the finite density states of interest. Another possibility is to separate the absolute value and the phase 
of the fermionic determinant and use the former to generate configurations and the latter in observables jfl) . 

At imaginary u. the measure remains positive and standard Monte Carlo techniques apply. The grand canonical 
partition function can be obtained by a Fourier transform 0, || . In this technique the dominant source of errors is the 
Fourier transform rather than the poor overlap. One can also use the fact that the partition function away from the 
transition line should be an analytic function of /Lt, and the fit for imaginary fi values could be analytically continued to 
real values of /i ||. At T sufficiently above the transition, both real and imaginary fi can be studied by dimensionally 
reducing QCD ]io| . Hamiltonian formulation may also help studying the problem jljj] . One can study adjoint matter 
and color superconductivity in two-color QCD |T^] . Nambu-Jona-Lasinio model was also used as an effective theory for 
strong interactions Jl3| . It has also been proposed to consider a for isospin rather than for baryon number || Q . 
Another approach avoids the sign problem by using cluster algorithms, in which negative against positive contributions 
are cancelled ]l5[ . 

We propose a method -an overlap improving multi-parameter reweighting technique- to reduce the overlap problem 
of the Glasgow method and determine the phase diagram in the T-fi plane. (Note, that a similar technique was successful 
for determining the phase diagram of the hot electroweak plasma jl6| e.g. on four-dimensional lattices, for which the 
applicability of a single-parameter reweighting was poor.) 

We study the system (say four-dimensional QCD with dynamical fermions) at Re(/i)=0 around its transition point. 
Using a Glasgow-type technique we calculate the determinants for each configuration for a set of /i, which, similarly to 
the Ferrenberg-Swendsen method can be used for reweighting. Using the average plaquette values of the individual 
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Figure 1: The average of the quark condensates at /3=5.085 as a function of Im(/i), for direct simulations (squares; their 
sizes give the errors), our technique (crosses) and Glasgow- type reweighting (dots). 



configurations for these partition functions an additional Ferrenberg-Swendsen ]17| reweighting can be done in the other 
parameter, thus in (3, too. For Re(^i)^0 and/or Im(/3)^0 significant cancellations of the complex phases appear, but 
exactly this feature is used in the determination of the zeros of the partition functions (Z), when looking for the transition 
point. Simultaneously reweighting in the two parameters /j, and we can keep the system on the transition line, which 
can be controlled e.g. by the inspection of the Lee- Yang zeros [[L8| of Z at complex (3. (Note, that the idea of performing 
a reweighting near the QCD critical line was already suggested in ||). This technique gives a good overlap with the 
original transition-like states. In principle any other parameter can be used for this type of reweighting. 

We test this method and illustrate its success compared to the Glasgow method. We present the exploratory results 
for the fj,-T phase diagram of the dynamical rtf—A staggered QCD. Simulations were done on L t — 4 lattices of 4 4 and 
4 • 6 3 . We estimate the phase diagram in physical units using the p mass (m p ) as the definition of the scale. 

Due to the small lattices and large spacings our estimate has systematic uncertainties, which can be reduced by 
approaching the continuum limit. The study of this limit is clearly not the goal of the present Letter. 

The Letter is organised as follows. Section II. presents the overlap ensuring multi-parameter reweighting technique. 
Section III. illustrates the applicability of the technique using rif—A dynamical QCD and gives the phase diagram in 
physical units. We conclude in Section IV. 



2 Overlap improving mult i- parameter reweighting 

Let us study a generic system of fermions ip and bosons (j), where the fermion Lagrange density is ipM^) 1 ^. Integrating 
over the Grassmann fields one gets: 

Z(a) = Jv(l)eM-Sboa(a,(f>)}detM(cj),a), (1) 

where a denotes a set of parameters of the Lagrangian. Eg. in the case of QCD with staggered quarks a consists of /3, 
the quark mass (m q ) and \i (which is included as exp(/i) and exp(— \x) multiplicative factors of the forward and backward 
timelike links, respectively). For some choice of the parameters a—a^ importance sampling can be carried out (e.g. for 
Re(/u)=0). Rewriting eq. (|l|) one obtains 

Z{a) = J V<j}cxp[-S bos (a , 4>)] det M(0, a ) 

(exp {-S bos (a, 0) + S bos (ao , 0)] ) • (2) 

L detM(</>, a ) J 

Now we treat the terms in the curly bracket as an observable -which is measured on each independent configuration- 
and the rest as the measure. It is well known that changing only one parameter of the set a the ensemble generated 
at ao provides an accurate value for some observables only for very high statistics. This is ensured by important but 
rare fluctuations as the mismatched measure occasionally sampled the regions where the integrand is large. This is the 
so-called overlap problem. Note however, that we have several parameters and the set a can be adjusted to ao to ensure 
much better overlap than obtained by varying only one, single parameter. Since the calculation of the determinants is 
expensive the most straightforward way to obtain a good overlap is to fix the parameters of the matrix M and adjust the 
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Figure 2: The average of the Polyakov lines (squares) and quark condensates (triangles) as a function of (3 at Im(/x) = 
0, Re(^i) = 0.3 (left panel). Histogram of the plaquettes at /3=4.938 and ^i=0.3 (right panel). The lattice volume is 4 4 . 



parameters which appear only in the bosonic action (in other words perform a Ferrenberg-Swendsen reweighting based 
on the bosonic part of the curly bracket). 

By simulating at a given set of bosonic couplings and using Ferrenberg-Swendsen reweighting we can get information 
on the system at other values of the couplings, even for complex ones. At finite volumes Z(a) has zeros -thus the free 
energy singularities- for complex values of these couplings. Standard finite size scaling techniques can be used to analyse 
the volume dependence of the Lee- Yang zeros. When looking for these zeros we use reweighting for some bosonic couplings. 
Simultaneously changing two -or more- parameters we can ensure that the system is reweighted along a transition line. 
This can be monitored by inspecting the Lee- Yang zeros of Z . 

3 Illustration: n/=4 dynamical QCD results 

For the case of lattice QCD at nonvanishing /i the above idea can be applied as follows. We will use two parameter 
reweighting, namely reweighting in (3 and /i. One performs the simulations at some /?, m and fi with Re(/i)=0 (note that 
purely imaginary fj, can be directly simulated). For each independent gauge configuration one calculates the average value 
of the plaquettes and the ratio of the determinants detM(//)/ det M(/i; Re(/i) = 0). For each // some 5(3 can be used 
to reweight with the measured plaquette values. By this way a much better overlap can be ensured than by reweighting 
only in /x (Glasgow method). 

We have tested these ideas in four-flavor QCD with m„=0.05 dynamical staggered quarks. The molecular dynamics 
Monte-Carlo code of the MILC collaboration was used |L{|. We checked that the determinants were calculated on 
independent configurations. Our statistical errors were obtained by a jackknife analysis. 

In order to directly check the applicability of our method we first collected 1200 independent V=4-6 3 configurations at 
Re(^)=lm(/x)=0 and used the Glasgow-reweighting and also our technique to study Re(/x)=0, Im(/Lt)^0. For our method 
we used the transition (3 (5.040) to generate the configurations, while for the Glasgow method /3=5.085 was used. At 
Re(/i)=0, lm(/x)^0 and at /3=5.085 direct simulations are possible. After performing these direct simulations as well, 
a clear comparison can be done. Figure |l| shows the predictions of the three methods for the quark condensate. The 
prediction of our method is in complete agreement with the direct results, whereas the predicition of the Glasgow-method 
is by several standard deviations off. Figure [l] indicates that the two-parameter reweighting is far more trustworthy than 
the single parameter one. Note, that imaginary chemical potential is a useful check on the proposed method, though it 
is different from a real chemical potential in biasing the ensemble towards non-zero baryon density. 

Based on these experiences we expect that our method is superior also at Re(^)^0. 

Next we study the physical Re(/i)^0 case. In order to have a better overlap and to check further the applicability of our 
reweighting technique we carried out simulations on 4 4 lattices at four different imaginary /i values lm(/i)=0, 0.1,0.2,0.3. 
The results obtained by the different runs are in complete agreement after reweighting. In the following we use our largest 
sample generated at lm(/i)=0. On this smaller V we used 13000 independent configurations. On the 4-6 3 lattice Im(/x)=0 
was used with 1200 independent configurations. The runs were carried out at the transition /3 values at Re(/z)=0. 

Let us illustrate that we are really at a transition point for nonvanishing fj,=0.3. Figure || shows the reweighted 
Polyakov-line and chiral condensate as a function of (3. Furthermore, we check the correctness of our Lee- Yang reweighting 
approach by showing the structure of a histogram. As it can be seen the turning point (indicating the coexistence of the 
two phases, thus the transition) is at » /3=4.94. At /i=0.3 the histogram method predicts a equal-weight double-peak 
stucture at /3 C =4. 939(5), which is in complete agreement wit the prediction of the Lee- Yang zero technique (see later). 
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Table 1: Lee- Yang zeros obtained at different /i values. 
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Figure 3: The phase diagram in the T-fi plane for nj — 4 QCD. The physical scale is set by m p . In physical units m q w 
25 MeV. The last point (/i «190 MeV) corresponds to our largest reweighted \i. 



Table 1. gives the real parts of the Lee- Yang zeros. Only zeros with the smallest imaginary parts are listed. The real 
parts of these zeros are usually used as a definition of the transition (3 at finite V. (Note, that the oo limit of the 
imaginary parts tells the order of the transition, cf. [|l8| ^|.) Based on V=4 4 and 4 ■ 6 3 we estimated the V—* oo limit 
by 1/V scaling. The critical (3 in this limit is used to transform the results to physical units. 

It is of particular interest to determine the phase diagram of QCD on the T-fi plane. Though the lattices we used 
are absurdly small and the spacing is quite large, it is illustrative to transform /3, m q and \x into physical units. Several 
parameters can be chosen to fix the scale, they give quite different values at our (3 couplings. In the present analysis we 
fixed the scale by m p =770 MeV. For small (3 values, studied by the present paper, m p can be obtained by interpolating 
between the strong coupling regime [^lf and the early measurements | j22| . Figure ^ shows the phase diagram in physical 
units. The errorbars indicate the statistical uncertainties reached on our sample of only O(10 3 — 10 4 ) configurations. 
Note, that L t —A lattices with the above definition of the scale restrict T to be larger than approximately 100 MeV (for 
clarity we used this value at the origin). 



4 Conclusions, outlook 

We proposed a method -an overlap improving multi-parameter reweighting technique- to numerically study non-zero fx 
and determine the phase diagram in the T-fi plane. We applied this technique for nf=A QCD with dynamical staggered 
quarks. We showed that for lm(/i)^0 the predictions of our method are in complete agreement with the direct simulations, 
whereas the Glasgow method suffers from the well-known overlap problem. Based on rather small statistics we were able 
to determine the critical gauge couplings as a function of the real chemical potential, which result was transformed into 
physical units. 

In this exploratory study we concentrated on the transition line separating the two phases. Clearly, the same technique 
can be applied to T-fi values somewhat below or above the line, for which the configurations should be collected at an 
appropriately chosen (3 -below or above the transition one- and at Re(/i)=0. 

Note, that the factor in curly braces in (2) is of order exp[— VS], with 6, of course, dependent on the parameters 
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a and and the configuration (f>. Making 6 small by a judicious choice of ag seems to work reasonably well at finite 
temperature. Nevertheless, there are two apparent limitations of the method. For large enough volumes the reweighting 
factor will always be exponentially suppressed, and at zero temperature 8 is larger than at finite T, thus the method most 
probably can not be applied to locate the transition at T=0. 

Our method can be easily applied to any number of Wilson fermions, whereas for H/=2 or nf =2+1 staggered fermions 
the situation is more complicated due to the ambiguity of the roots of the determinants p0| . 

Note, that the present reweighting does not provide a general solution to the sign problem in QCD, but could be 
useful to locate the critical point of QCD, for which we are interested in a relatively small /x and high T region of the 
phase diagram and may get away with relatively small volumes. Another application is to determine the curvature of the 
phase diagram p3| ]. 

We thank F. Csikor and I. Montvay for useful discussions and comments on the manuscript. This work was partially 
supported by Hungarian Science Foundation grants No. OTKA-T37615/T34980/T29803/T22929/M28413-37071/OM- 
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